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The transition to turbulence exhibits remarkable spatio-temporal behavior that continues to defy 
detailed understanding. Near the onset to turbulence in pipes, transient turbulent regions decay 
either directly or, at higher Reynolds numbers through splitting, with characteristic time-scales 
that exhibit a super-exponential dependence on Reynolds number. Here we report numerical sim¬ 
ulations of transitional pipe flow, showing that a zonal flow emerges at large scales, activated by 
anisotropic turbulent fluctuations; in turn, the zonal flow suppresses the small-scale turbulence 
leading to stochastic predator-prey dynamics. We show that this “ecological” model of transi¬ 
tional turbulence reproduces the super-exponential lifetime statistics and phenomenology of pipe 
flow experiments. Our work demonstrates that a fluid on the edge of turbulence is mathematically 
analogous to an ecosystem on the edge of extinction, and provides an unbroken link between the 
equations of fluid dynamics and the directed percolation universality class. 


Introduction. Fluids in motion are generally found in 
one of two generic states. The most common - turbu¬ 
lence is found at sufficiently large characteristic speeds 
U, depending on the kinematic viscosity v and the char¬ 
acteristic system scale D ; turbulent flows are complex, 
stochastic, and unpredictable in detail. At lower veloci¬ 
ties, the fluid is said to be laminar: its flow is simple, de¬ 
terministic and predictable. In between these two states, 
conventionally delineated by the dimensionless control 
parameter known as the Reynolds number Re = UD/v , 
is a transitional regime that occurs for 1700 < Re < 2300 
in pipes, and which has presented a challenge to experi¬ 
ment and theory since Osborne Reynolds’ original obser¬ 
vation of intermittent “flashes” of turbulence (T. Today, 
Reynolds’ flashes are known as puffs [Tj, and their be¬ 
havior has been characterized very precisely through a 
series of physical and numerical experiments performed 
during the last decade or so [3HS] (for a recent review, 
see H) culminating in the tour de force observation of 
a super-exponential functional dependence of the life¬ 
time T d of puffs as a function of Re [8] : lnlnr d oc Re. 
For Reynolds numbers based upon pipe diameter D of 
around 2300 , turbulence is sustained longer than the abil¬ 
ity to observe its lifetime in finite systems, and the puffs 
become unstable through a new dynamical processes in 
which the leading edge breaks away and nucleates the for¬ 
mation of a new puff some distance downstream [HHia. 
The puff-splitting occurs on a characteristic time that 
decays super-exponentially with increasing Re. These 
phenomena are likely to be generic. For example, super¬ 
exponential scaling behavior near the transition to tur¬ 
bulence has also been reported in plane Poiseuille flow 
[ 13 ] and Taylor-Couette flow [ 14 ] . In addition, optical 
analogues of the laminar-turbulence transition have been 
observed in laser-driven optical fibers EM- 

The theoretical account of these phenomena has fo¬ 
cused primarily on the existence and interactions between 
nonlinear solutions of the Navier-Stokes equations, peri¬ 
odic orbits and streamwise vortices [T6U2T] . and the dy¬ 
namics of long-lived chaotic transients [22] . An alterna¬ 


tive line of inquiry has been to characterize the statisti¬ 
cal properties of transitional turbulence through ad hoc 
model equations. These have been motivated either by 
perceptive analogies to excitable media [23] or by phase 
transition universality arguments that begin with the no¬ 
tion that the laminar state is an absorbing one [24], and 
show quantitatively how super-exponential decay results 
from the generic universality class for non-equilibrium 
absorbing processes [23]: directed percolation (DP) [23] , 
Both approaches reflect an important aspect of the dy¬ 
namics, namely that a certain minimum level of energy 
is required to sustain turbulent puffs [27], leading to a 
further connection with extreme value statistics [28]. 

It is the statistical behavior near the transition which 
concerns us here: how do the various spatial-temporal 
modes that are excited give rise to such remarkable life¬ 
time statistics? What is the universality class of this 
transition, in terms of its fluctuation characteristics? Are 
there simplified effective descriptions that bridge the gap 
between the underlying fluid dynamics and the large- 
scale statistical properties? And how do these emerge 
from the underlying Navier-Stokes equations that govern 
all hydrodynamic phenomena? 

The purpose of this report is to address these questions 
using the same approach by which phase transitions are 
understood in condensed matter physics [29]. There it 
is well-established that universal aspects of phase tran¬ 
sitions, such as the phase diagram, critical exponents, 
scaling functions etc. are all described quantitatively 
by an effective coarse-grained theory (“Landau theory”) 
that contains the symmetry-allowed collective and long- 
wavelength modes, without requiring excessive realism 
at the microscopic level of description. Being based on 
symmetry principles, the individual symmetry-allowed 
terms in Landau theory do not require detailed deriva¬ 
tion from the microscopic level of description. This is 
fortunate given that there is usually no good, uniformly 
valid approximation scheme to derive formally and sys¬ 
tematically these terms and their coefficients from first 
principles. 
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For the transitional turbulence problem, the analyti¬ 
cal difficulties are acute. Therefore, in this paper, we 
avoid approximations which are difficult to justify sys¬ 
tematically. Instead, we use direct numerical simulation 
to identify the important collective modes which exhibit 
an interplay between large-scale fluctuations and small- 
scale dynamics at the onset of turbulence, and thence to 
write down the corresponding minimal stochastic model, 
in the spirit of the Landau theory of phase transitions. 
As we will show, this coarse-grained description is able to 
account for the principal experimental observations and 
to predict quantitatively the puff lifetime and splitting 
behavior observed near the transition. Our stochastic 
model can be transformed using standard statistical me¬ 
chanics techniques into a field theory known to be in 
the universality class of directed percolation. These re¬ 
sults provide an unbroken link between the equations of 
fluid dynamics and the directed percolation universality 
class JH (2B|. Our approach is thus a precise parallel 
to the way in which phase transitions are understood in 
condensed matter physics, and shows that concepts of 
universality and effective theories are applicable to the 
laminar-turbulence transition. 

Observation of predator-prey dynamics in Navier- 
Stokes equations. To address these questions without 
making questionable or non-systematic approximations, 
we have performed direct numerical simulations (DNS) 
of the Navier-Stokes equations in a pipe, using the open- 
source code “Open Pipe Flow” [3D], as described in Ap¬ 
pendix [A] We denote the time-dependent velocity devia¬ 
tion from the Hagen-Poiseuille flow by u = (u z ,ug,u r ). 
Because we were interested in transitional behavior, we 
looked for a decomposition unnuEiEa of large-scale 
modes that would indicate some form of collective be¬ 
havior, as well as small-scale modes that would be repre¬ 
sentative of turbulent dynamics In particular, we report 
here the behavior of the velocity field (u z ,ug,u r ), where 
the bar denotes average over z and 9 , and u r = 0. We 
refer to this as the zonal flow (ZF). In Fourier space, the 
zonal flow is given by u{k = 0, m = 0, r), where k is the 
axial wavenumber and m is the azimuthal wavenumber, 
r is the real space radial coordinate and the tilde denotes 
Fourier transform in the 9 and z directions only. Turbu¬ 
lence was represented by short-wavelength modes, whose 
energy is E T {t) = § E| fc |>i, H >i / I u(k, m, r) | 2 dF. 

Shown in Figure 0 A ) is a time series for the energy 
EzFit) of the zonal flow, compared with the energy Exit) 
of the turbulent energy. The curves show clear persis¬ 
tent oscillatory behavior, modulated by long-wavelength 
stochasticity as shown in the phase portrait of Figure 
0B)- In Figure QC), we have calculated the phase shift 
between the turbulence and zonal flows, with the result 
that the turbulent energy leads the zonal flow energy by 
~ 7t/2. This suggests that these oscillations can be inter¬ 
preted as a time-series resulting from activator-inhibitor 
dynamics, such as occurs in a predator-prey ecosystem. 
Predator-prey ecosystems are characterized by the fol¬ 
lowing behavior: the “prey” mode activates the “preda- 
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FIG. 1. Predator-prey oscillations in transitional turbu¬ 
lent pipe flow. (A) Energy vs. time for the zonal flow (or¬ 
ange) and turbulent modes (green). (B) Phase portrait of 
the zonal flow and turbulent modes as a function of time, 
with color indicating the earliest time in dark blue progress¬ 
ing to the latest time in light green. (C) Phase shift be¬ 
tween the turbulent and zonal flow modes as a function of 
frequency, showing that the turbulence leads the zonal flow 
by 7r/2 consistent with predator-prey dynamics. The phase 

shift 9 (to) = tan -1 (^Im[C(to)]/Re[C(oj)]^ and is shifted to 

be positive, where C(u>) is the Fourier transform of the corre¬ 
lation function between the turbulence and the zonal flow in 
(A). The red line corresponds to the dominant frequency in 
the power spectrum. The phase shift near small oj is scatter 
due to the finite time duration of the time series. (D) Veloc¬ 
ity field configuration of the zonal flow mode u. The color 
bar indicates the value of u z . (E) Snapshot of the Reynolds 
stress gradient and zonal flow time-derivative as functions of 
r. (F) Reynolds stress gradient and zonal flow time derivative 
as functions of time. The evident proportionality shows that 
zonal flow dynamics is driven by the radial gradient of the 
Reynolds stress. 
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tor” mode, which then grows in abundance. At the same 
time, the growing predator mode begins to inhibit the 
prey mode. The inhibition of the prey mode starves the 
predator mode, and it too becomes inhibited. The in¬ 
hibition of the predator mode allows the prey mode to 
re-activate, and the population cycle begins again. 

The flow configuration for the predator mode is shown 
in Figure [T] (D), and consists of a series of azimuthally 
symmetric modes with direction reversals as a function 
of radius r. Such banded shear flows are known as zonal 
flows and are of special significance in plasma physics, as- 
trophysical and geophysical flows, owing to their role in 
regulating turbulence |33j . The purely azimuthal com¬ 
ponent of the zonal flow, denoted by ug(r) is spatially 
uniform in z , and the lack of a radial component means 
that it is not driven by pressure gradients. Thus it can 
only exist due to nonlinear interactions with turbulent 
modes. In this sense, it is a collective mode, one with 
special significance for transitional turbulence. 

The simplest way for such an azimuthal shear flow to 
couple to turbulent fluctuations is through the Reynolds 
stress t: however, a uniform Reynolds stress cannot drive 
a shear flow, so the first symmetry-allowed possibility is 
the radial gradient of the Reynolds stress 33], as ex¬ 
pressed in the Reynolds momentum equation. Thus, to 
probe the dynamics that govern the emergence of the 
zonal flow, we have calculated the time-averaged radial 
gradient of the instantaneous Reynolds stress, r = u' g -u' r , 
where u'(z, 9,r ) = u — u, and show in Figure [l] (F) the 
4.5-time-unit-running-mean time series of —dtue and the 
radial gradient d r r. Both quantities have been averaged 
over 0 < z < L, 0 < 9 < 2n and i? 0 < r < R, where 
R 0 = 0.6411?, and the resulting time series are clearly 
highly correlated. 

In general, it is the case that zonal flows are driven by 
statistical anisotropy in turbulence, but are themselves 
an isotropizing influence on the turbulence through their 
coupling to the Reynolds stress [34H36] , The fact that 
turbulence anisotropy activates the zonal flow, and that 
zonal flow inhibits the turbulence is responsible for the 
predator-prey oscillations observed in the numerical sim¬ 
ulations. 

Lifetime of stochastic predator-prey populations. 
Phase transition theory [29] would suggest that the emer¬ 
gence of a zonal flow collective mode dominates the non¬ 
equilibrium transition of pipe flow from the laminar to 
the turbulent state, through the predator-prey interac¬ 
tion with the small scale velocity fluctuations. Such 
a “two fluid” effective field description of transitional 
turbulence implies that stochastic predator-prey popula¬ 
tions should undergo spatio-temporal fluctuations whose 
functional form matches precisely the observations for 
the lifetime and splitting time of turbulent puffs in a 
pipe. To test this idea, we have performed simulations 
of a spatially-extended stochastic predator-prey ecosys¬ 
tem, in a quasi-one-dimensional geometry to mimic the 
pipe environment. The specific system has three trophic 
levels: nutrient (E), Prey (B) and Predator (A), which 


correspond in the fluid system to laminar flow, turbulence 
and zonal flow respectively. The interactions between in¬ 
dividual representatives of these levels are given by the 
following rate equations 

Ai Ei, Bi^ E it Ai + B, A, + A,, 

(ij) 

Bi + Ej Bi + B h Bi A Ai, 

(ij) 

Ai + Ej -> Ej + Aj , Bj + Ej -> Ej T Bj. (1) 

(ij) (ij) 

where 9a and ds are the death rates of A and B, p is the 
predation rate, b is the prey birth rate due to consump¬ 
tion of nutrient, (ij) denotes hopping to nearest neighbor 
sites, D is the nearest-neighbor hopping rate, and to is 
the point mutation rate from prey to predator, which 
models the induction of the zonal flow from the turbu¬ 
lence degrees of freedom. 

We are primarily interested in long-wavelength prop¬ 
erties of the system, at least in the vicinity of the 
turbulence transition, where we expect the transverse 
correlation length to be larger than the pipe diame¬ 
ter, implying that the behavior is in fact quasi-one- 
dimensional. The crossover phenomena associated with 
this have been discussed previously [26], and thus our 
quasi-one-dimensional model should be appropriate and 
quantitatively correct near the transition. 

In our simulation, described in Appendix |Bj the con¬ 
trol parameter is the prey birth rate b. When b is small 
enough, the population is metastable, and cannot sustain 
itself: all individuals, both predator and prey, eventually 
die within a finite lifetime T d (b). As b increases, the life¬ 
time of the population increases rapidly: in particular the 
prey lifetime increases rapidly with b. At large enough 
values of b, the decay of the initial population is not ob¬ 
served, but instead the initially localized population pro¬ 
liferates, spreading outwards and spontaneously splitting 
into multiple clusters, as shown in the space-time plot of 
clusters of prey of Figure [2] (A). 

To quantify these observations, we have measured both 
the lifetime of population clusters in the metastable re¬ 
gion and their splitting time using a procedure directly 
following that of the turbulence experiments and simu¬ 
lations HD, and described in Appendix [C] We comment 
that both timescales involve implicitly measurements of 
quantities that exceed a given threshold, and thus it is 
natural that the results are found to conform to extreme 
value statistics [26] GHj ■ 

In Figure [2] (A) we show the phenomenology of the 
dynamics of initial clusters of prey, corresponding to the 
predator-prey analogue for the experiments in pipe flow 
which followed the dynamics of an initial puff of turbu¬ 
lence injected into the flow [5]. Depending upon the prey 
birth rate, the cluster decays either homogeneously or by 
splitting, precisely mimicking the behavior of turbulent 
puffs as a function of Reynolds number. The extraction 
from data of decay times is described in Appendix [C] In 
Figure [2] (B) is shown the semi-log plot of lifetime for 






4 




0.02 0.03 0.04 0.05 0.06 0.07 


b 


FIG. 2. Stochastic predator-prey model reproduces the phe¬ 
nomenology of transitional pipe turbulence. Decay lifetime, 
t , and splitting time, t s , of clusters of prey are memory¬ 
less processes and obey super-exponential statistics as a func¬ 
tion of prey birth rate. To compare with the experiments 
[8], predator-prey dynamics are performed in two-dimensional 
pipe geometry as described in the text. (A) World line of clus¬ 
ters of prey splitting to form predator-prey traveling waves. 
The color measures the local density of prey, corresponding 
to intensity of turbulence in pipe flow. In the simulation, 
the dimensionless parameters are D = 0.1, b = 0.1, p = 0.2, 
d,A = 0.01, d B = 0.01 and m = 0.001. In the model simu¬ 
lated, diffusion is isotropic, not biased as would be the case 
corresponding to a mean flow, where the clusters will accumu¬ 
late at large times with a well-defined separation set by the 
depletion zone of nutrient behind each predator-prey travel¬ 
ing wave. (B) Log decay lifetime of prey cluster and splitting 
time as a function of prey birth rate. The upward curva¬ 
ture signifies super-exponential behavior. The parameters are 
D = 0.01, p = 0.1, d A = 0.015, d B = 0.025 and m = 0.001. 
Inset: Double log lifetime vs prey birth rate, showing the fit 
to the following functional forms: the dashed curve is given 
by T d /ro = exp(exp(46.5396 — 0.731)), and the solid curve is 
given by t 3 /tq = exp(exp(—31.1486 — 3.141)). 


both decay and splitting as a function of prey birth rate, 
the upward curvature indicative of super-exponential be¬ 
havior. The inset to Figure [2] (B) shows a double expo¬ 
nential plot of puff lifetime and splitting time vs. prey 
birth rate, the straight line being the fit to the functional 
form indicated in the caption. These figures indicate a 


remarkable similarity to the corresponding plots obtained 
for transitional pipe turbulence in both experiments [8] 
and direct numerical simulations m, and demonstrate 
conclusively that experimental observations are well cap¬ 
tured by an effective two-fluid model of pipe flow turbu¬ 
lence with predator-prey interactions between the zonal 
flow and the small scale turbulence. 

Universality class of the laminar-turbulence transi¬ 
tion in pipes. The two-fluid predator-prey model ex¬ 
pressed by Equations 0 exhibits a rich phase dia¬ 
gram that captures the main features observed in tran¬ 
sitional turbulence in pipes. The transition to puff¬ 
splitting can be identified with a change of stability 
of the spatially-uniform mean-field predator-prey coex¬ 
istence point, where a stable node becomes a stable fo¬ 
cus or spiral with increasing birth rate. In the language 
of predator-prey systems, this corresponds to the break¬ 
down of spatially homogeneous prey domains into peri¬ 
odic traveling wave states. The phase diagram is sketched 
in Figure [3] along with the corresponding phase diagram 
for transitional pipe turbulence as determined by exper¬ 
iment. The phenomenology of the predator-prey system 
mirrors that of turbulent pipe flow. 

In order to determine the universality class of the non¬ 
equilibrium phase transition from laminar to turbulent 
flow, we use the two-fluid predator-prey mode in Equa¬ 
tions |l]) . Near the transition to prey extinction, the prey 
population is very small and no predator can survive, and 
thus Equations 0 simplify to 

Bi + Ej ——>■ Bi + Bj, 

(U) 

Bi + Ej - > Ei + Bj . ( 2 ) 

(*?> 

These equations are exactly those of the reaction- 
diffusion model for directed percolation [37] . A more de¬ 
tailed and systematic way to reach this conclusion is to 
represent Equations [l] exactly in path integral form using 
the Doi formalism [5TH42] . The resulting action simpli¬ 
fies near the transition to that of Reggeon field theory 
S3®, which has been shown to be in the universality 
class of directed percolation [25105]. Numerical simu¬ 
lations of 3 + 1 dimensional directed percolation in a 
pipe geometry have reproduced the statistics and behav¬ 
ior of turbulent puffs and slugs in pipe flow [TQ 0B] , and 
a detailed comparison between theory and experiment 
has been presented 02 ■ The super-exponential behavior 
of DP might seem to contradict the expectation based 
upon the known critical behavior (e.g., see Ref. [48]). 
However, it is important to recognize that the usual ex¬ 
ponents relate to DP starting from a single seed, whereas 
the experiments and simulations are conducted with an 
extended seed that has a finite length or number of seed 
points. These points behave as independent identically- 
distributed random variables as long as the transverse 
correlation length is much smaller than the seed size, but 
once the correlation length is of order the seed size, the 
seed is effectively a single correlated extended source, and 
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FIG. 3. Schematic phase diagram for transitional pipe tur¬ 
bulence as a function of Reynolds number compared with the 
phase diagram for predator-prey dynamics as a function of 
prey birth rate. Above each phase is shown a typical flow or 
predator-prey configuration, indicating the similarity between 
the turbulent pipe and ecosystem dynamics. 


once the correlation length is much larger than this size, 
there will be a crossover to the usual DP exponents. 

Discussion. The observation of the emergence of a 
zonal flow, excited by the developing turbulent degrees of 
freedom and the demonstration of its role in determining 
the phenomenology of transitional pipe turbulence has an 
interesting consequence: the zonal flow can be assisted by 
rotating the pipe, and this should catalyze the transition 
to turbulence, causing it to occur at lower Re. Indeed 
experiments on axially-rotating pipes [49] are consistent 
with this prediction. 

The idea that predator-prey dynamics can arise in tur¬ 
bulence is by no means new, and such behavior was pro¬ 
posed by Diamond and collaborators [331 EES ;5TJ many 
years ago in the context of the interaction between drift- 
wave turbulence and zonal flows in tokomaks; indeed the 
predator-prey oscillations were recently observed in toko¬ 
maks [52H5B1 and in a table-top electroconvection ana¬ 
logue of the L-H transition {35]. The new ingredient we 
have presented in this paper is the observation of a zonal 
flow emerging in the very simple setting of transitional 
pipe turbulence, and the demonstration of predator-prey 
oscillations that suggest a minimal mesoscale model of 
transitional turbulence that accounts for the observed 
statistical behavior of puffs. The logic of our work is that 
we used DNS to identify the important collective modes 
at the onset of turbulence—the predator-prey modes— 
and then wrote down the simplest minimal stochastic 
model to account for these observations. This model 
predicts without using the Navier-Stokes equations the 


puff lifetime and splitting behavior observed in experi¬ 
ment. This approach is a precise parallel to that used 
in the conventional theory of phase transitions, where 
one builds a Landau theory, a coarse-grained (or effec¬ 
tive) theory, using symmetry principles. This intermedi¬ 
ate level description can then be used as a starting point 
for renormalization group analysis to compute the critical 
behavior. 

Note that it is neither necessary nor desirable to derive 
the coarse-grained effective theory from the microscopic 
description. Such a derivation is not necessary, because 
any such derivation would need a small parameter and 
would thus only have limited validity due to the analyt¬ 
ical approximations made. Instead, the effective theory 
can be obtained from symmetry principles. A familiar 
example of this situation is that even though the Navier- 
Stokes equations can be derived from Boltzmann’s ki¬ 
netic equations for gases, such a derivation would imply 
that the Navier-Stokes equation is only valid for dilute 
gases. In fact, the Navier-Stokes equations are an ex¬ 
cellent description for dense liquids as well, and can be 
obtained by perfectly satisfactory phenomenological and 
symmetry arguments. The derivation from Boltzmann’s 
kinetic theory is inherently limited by the regime of va¬ 
lidity of the kinetic theory— low density—and this leads 
to an unnecessarily restrictive derivation of the equa¬ 
tions of fluid dynamics. The reason why an analytical 
derivation of the coarse-grained theory is unnecessary is 
that even if the coefficients of the terms could be com¬ 
puted in the order parameter expansion of the Landau 
theory, they do not come into the exponents or scal¬ 
ing functions anyway, and thus they do not affect the 
critical behavior. In the present context, the effective 
theory we obtain is then mapped exactly into directed 
percolation. The observation of an emergent zonal flow 
and predator-prey oscillations with its attendant minimal 
“two-fluid” model, provide a direct and unbroken link 
between the Navier-Stokes equations and the directed 
percolation universality class for transitional turbulence. 
Our work underscores not only the potential importance 
of zonal flows in other transitional turbulence situations, 
but also shows the utility of coarse-grained effective mod¬ 
els for non-equilibrium phase transitions, even to states 
as perplexing as fluid turbulence. 
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SUPPLEMENTARY MATERIAL 


Appendix A: Direct numerical simulations of the 
Navier-Stokes equations. 

We performed direct numerical simulations (DNS) of 
the Navier-Stokes equations in a pipe, using the open- 
source code “Open Pipe Flow” [30]. The equations were 
solved using a pseudo-spectral method in cylindrical co¬ 
ordinates |30j . having 60 grid points in the radial (r) di¬ 
rection, 32 Fourier modes in the azimuthal (9) direction 
and 128 modes in the axial (z) direction. Such a model 
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is of course a reduced description of reality, but the main 
features can be well captured, with a slight renormal¬ 
ization of the Re needed to compare with experiment 
m■ For the accuracy required in our study, we used 32 
modes in the azimuthal direction, typical of most stud¬ 
ies in the literature, although we note that useful, semi- 
quantitative results have been obtained with as few as 
2 modes [30,. The spatial resolutions were chosen such 
that the resolvable power spectra span over six orders of 
magnitude. The pipe length L is 20 times its radius i?, 
with periodic boundary conditions in the z direction [30]. 
With this resolution, the transition to turbulence occurs 
in a range of Re numbers between 2200 and 3000, and 
moves to smaller Re at still higher resolution. We re¬ 
port here measurements at Re = 2600, slightly above the 
transition m- The mass flux and Re = 2600 were held 
constant in time [30] . The laminar flow is the Hagen- 
Poiseuille flow, which was independent of time as the 
mass flux was held constant [30] . 


Appendix B: Stochastic simulations of predator-prey 
dynamics. 



. Q . ... . 
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The specific system has three trophic levels: nutrient 
(E), Prey (B) and Predator (A), which correspond in 
the fluid system to laminar flow, turbulence and zonal 
flow respectively. Such a system can be naively modeled 
by the Lotka-Volterra ordinary differential equations 153- 
155] , which in the case of ecosystems with finite resources 
do not permit long-time persistent oscillatory solutions, 
unless additional biological details such as functional re¬ 
sponse are included. In fact, it is necessary to include 
the dynamics of individual birth-death events, and when 
this is done correctly, it is found that the number fluc¬ 
tuations drive the population oscillations m through 
resonant amplification. Thus, we use a stochastic model 
at the outset. 

The interactions between individual representatives 
of these levels are given by Equations [l] We simulated 
these equations on a 401 x 11 lattice in two dimensions, 
intended to emulate the pipe geometry. Lattice sites 
were only allowed to be occupied by one of E, A or B. 
The predator (A) and prey (B) are additionally allowed 
to diffuse via random walk on the lattice with diffusion 
coefficient 0.1 in units of the square of the lattice spacing 
divided by the time step (set equal to unity). The initial 
conditions for the simulations were a random popula¬ 
tion of prey and predator, occupying with probability 
4/5 and 1/5 respectively on the lattice sites between 
x £ [—15,15] and y £ [—5, 5] where x labels the direction 
along the axis of the ecosystem (pipe) and y labels the 
transverse direction. The predator-prey dynamics in Eq. 
[T] was implemented by the following algorithm: at each 
time step, a site i is randomly chosen, a random number 
s, is generated from the uniform distribution between 
zero and one. The behavior on the site is decided by 
the random number: (1) if s < 1/6 and the site i is 


FIG. SI. Stochastic predator-prey model reproduces the 
phenomenology of transitional pipe turbulence. Decay life¬ 
time, t , and splitting time, t s , of clusters of prey are mem¬ 
oryless processes and obey super-exponential statistics as a 
function of prey birth rate. To compare with the experiments 
[8], predator-prey dynamics are performed in two-dimensional 
pipe geometry as described in the text. The dimensionless pa¬ 
rameters in the simulation are D = 0.01, p = 0.1, d,A = 0.015, 
ds = 0.025 and m = 0.001. (A) Log survival probability of 
prey cluster vs. time during homogeneous decay to extinc¬ 
tion. Here the characteristic time scale that is estimated by 
ro ~ 200. (B) Log survival probability of prey cluster vs. time 
during decay to splitting. (C) Survival probability of prey 
cluster as a function of prey birth rate during homogeneous 
decay to extinction. (D) Survival probability of prey cluster as 
a function of prey birth rate during decay to splitting. (E) Log 
inverse lifetime of prey cluster, as a function of prey birth rate 
during homogeneous decay to extinction (left curve, r d ) and 
during decay to splitting (right curve, r a ). The dashed curve 
is given by ro/r d = 1/ exp(exp(46.539b—0.731)), and the solid 
curve is given by tq/t s = 1/ exp(exp(—31.148b — 3.141)). 


occupied by any individual, and if a randomly chosen 
neighbor site is empty, then that individual diffuses to 
the random neighboring site with rate y = 0.01 (i.e. this 
reaction happens if another uniformly distributed ran¬ 
dom number is less than 1 — exp /x); (2) if 1/6 < s < 1/3 
and the site i is occupied by a prey individual, and if 
a randomly chosen neighbor site, j, is empty, then one 
prey individual is born on the site j with rate b; (3) if 
1/3 < s < 1/2 and the site i is occupied by a predator 
individual, and if a randomly chosen neighbor site, j, is 
occupied by a prey individual, then the prey individual 
is replaced by a new-born predator individual with rate 













































p; (4) if 1/2 < s < 2/3 and the site i is occupied by a 
predator individual, that predator individual dies with 
rate (5) if 2/3 < s < 5/6 and the site i is occupied 
by a prey individual, that prey individual dies with rate 
ds\ (6) if 5/6 < s < 1 and the site i is occupied by a 
prey individual, then the prey individual is replaced by a 
predator individual with rate m. Then within the same 
time step, the above processes are repeated 401 x 11 
times so that on average one reaction takes place at each 
lattice site in the system. 


Appendix C: Measurement of decay and splitting 
lifetimes. 

We measured both the lifetime of population clusters 
in the metastable region and their splitting time using 
a procedure directly following that of the turbulence ex¬ 
periments and simulations m- To this end, we moni¬ 
tor the coarse-grained prey population density ns(i) = 
+ 3,1)/{H + 1)/(2J + 1) - 0.25, 
where H is the height of the system (11 lattice units) 
and J = 3. The lifetime of prey clusters is defined as 
the time it takes for the last prey individual to die. The 
cluster splitting time is defined as the first time that the 
distance between the edges of two coarse-grained prey 
clusters exceed 25 unit sites. We comment that both 


timescales involve implicitly measurements of quantities 
that exceed a given threshold, and thus it is natural that 
the results are found to conform to extreme value statis¬ 
tics [2BJ [55]. 


In Figure SI we show the phenomenology of the dy¬ 
namics of initial clusters of prey, corresponding to the 
predator-prey analogue for the experiments in pipe flow 
which followed the dynamics of an initial puff of turbu¬ 
lence injected into the flow [8]. Depending upon the prey 
birth rate, the cluster decays either homogeneously or 
by splitting, precisely mimicking the behavior of turbu¬ 


lent puffs as a function of Reynolds number. Figure SI 


(A) and (B) show that the decay is exponential in time, 
indicating that it is a memoryless process with a single 
time constant. Figure [ST] (C) and (D) show that the 
survival probability is a sigmoidal curve, whose inverse 
lifetime as a function of prey birth rate is plotted in a log 
scale in Figure [ST](E). If the lifetime were an exponential 
function, this curve would be a straight line with nega¬ 
tive slope. The downward curvature is a manifestation of 
super-exponential behavior. These figures indicate a re¬ 
markable similarity to the corresponding plots obtained 
for transitional pipe turbulence in both experiments [5] 
and direct numerical simulations EH, and demonstrate 
conclusively that experimental observations are well cap¬ 
tured by an effective two-fluid model of pipe flow turbu¬ 
lence with predator-prey interactions between the zonal 
flow and the small scale turbulence. 




